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We propose a systematic procedure for constructing effective lattice fermion models for narrow- 
band compounds on the basis of first-principles electronic-structure calculations. The method is 
illustrated for the series of transition-metal (TM) oxides: SrV03, YTi03, V2O3, and Y2M02O7, 
whose low-energy properties are linked exclusively to the electronic structure of an isolated tig 
band. The method consists of three parts, starting from the electronic structure in the local-density 
approximation (LDA). (i) construction of the kinetic-energy Hamiltonian using formal downfolding 
method. It allows to describe the band structure close to the Fermi level in terms of a limited number 
of (unknown yet) Wannier functions (WFs), and eliminate the rest of the basis states, (ii) solution 
of an inverse problem and construction of WFs for the given kinetic-energy Hamiltonian. Here, we 
closely follow the construction of the basis functions in the liner-mufnn-tin-orbital (LMTO) method, 
and enforce the orthogonality of WFs to other band. In this approach, one can easily control the 
contributions of the kinetic energy to the WFs. (iii) calculation of screened Coulomb interactions 
in the basis of auxiliary WFs. The latter are defined as the WFs for which the kinetic-energy 
term is set to be zero. Meanwhile, the hybridization between TM d and other atomic states is well 
preserved by the orthogonality condition to other bands. The use of auxiliary WFs is necessary 
in order to avoid the double counting of the kinetic-energy term, which is included explicitly in 
the model Hamiltonian. In order to calculate the screened Coulomb interactions we employed a 
hybrid approach. First, we evaluate the screening caused by the change of occupation numbers 
and the relaxation of the LMTO basis functions, using the conventional constraint-LDA approach, 
where all matrix elements of hybridization connecting the TM d orbitals and other orbitals are set 
to be zero. Then, we switch on the hybridization and evaluate the screening of on-site Coulomb 
interactions associated with the change of this hybridization in the random-phase approximation. 
The second channel of screening appears to be very important, and results in relatively small value of 
the effective Coulomb interaction for isolated t2 g bands (about 2-3 eV, depending on the material) . 
We discuss details of this screening and consider its band-filling dependence, frequency dependence, 
influence of the lattice distortion, proximity of other bands, as well as the effect of dimensionality 
of the model Hamiltonian. 

PACS numbers: 71.10.Fd; 71. 15. Mb; 71.28.+d; 71.15.Ap 



I. INTRODUCTION 

Many successes of modern solid-state physics and chemistry are related with the development of the Hohenberg- 
Kohn-Sham density- functional theory (DFT), 1,2 which is designed for the ground state and based on the minimization 
of the total energy E[p) with respect to the electron density p. For practical applications, DFT resorts to iterative 
solution of single-particle Kohn-Sham (KS) equations, 

-^y 2 + V H + V xc + V cxt )i Jl =e^ l , (1) 



2m 

together with the equation for the electron density: 

p = 5>^| 2 , (2) 

i 

defined in terms of eigenfunctions (ipi), eigenvalues (e^), and the occupation numbers (n^) of KS quasipar- 
ticles. Different terms in Eq.QJ are correspondingly the kinetic-energy operator, the Hartree potential, the 
exchange-correlation potential, and the external potential. In the following we will also reserve the notation 
Hks = — (fr 2 /2?™)V 2 +VH+Vxc+Kxt f° r tne total KS Hamiltonian in the real (r) space. 

The exchange-correlation potential is typically treated in the local-density approximation (LDA). It employs an 
analytical expression for Vxc [p] borrowed from the theory of homogeneous electron gas in which the density of the 
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electron gas is replaced by the local density of the real system. LDA is far from being perfect and there are many 
examples of so-called strongly-correlated materials for which the conventional LDA appears to be insufficient both for 
the excited-state and ground-state propertiesi ?' 4 '^ 

A typical situation realized in transition-metal (TM) oxides is shown in Fig. ^ We would like to emphasize two 
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FIG. 1: Total and partial densities of states for S1-VO3, YT1O3, V2O3, and Y2M02O7 in the local-density approximation. The 
shaded area shows the contributions of transition-metal d-states. Other symbols show positions of the main bands. The Fermi 
level is at zero energy. 

points. 

(i) The common feature of many TM oxides is the existence of the well isolated narrow band (or the group of bands) 
located near the Fermi level and well isolated from the rest of electronic stares. For compounds shown in Fig. ^ this 
is the TM ti g band, which is sandwiched between 0(2p) band (from below) and a group of bands (from above), which 
have an appreciable weight of the TM e g states (the meaning of many notations will become clear in Sec. IVII where 
we will discuss details of the crystal and electronic structure for the considered oxide compounds). Electronic and 
magnetic properties of these compounds are predetermined mainly by the behavior of this ti g band. The effect of 
other bands can be included indirectly, through the renormalization of interaction parameters in the ti g band. 

(ii) The LDA description appears to be especially bad for the ti g states located near the Fermi level. It often fails to 
reproduce the insulating behavior of these compounds, as well as the correct magnetic ground state, which is directly 
related with the existence of the band gap£ The source of the problem is know to be the on-site Coulomb correlations, 
whose form is greatly oversimplified in the model of homogeneous electron gas. Therefore, the basic strategy which 
was intensively pursued already for more than decade was to incorporate the physics of on-site Coulomb correlations in 
LDA and to solve this problem using modern many-body techniques. This way of thinking gave rise to such directions 
as LDA+U (e.g., Refs.HH and0) and LDA+DMFT (dynamical mean-field theory, Ref.||). 

Taking into account two above arguments, we believe that the most logical way to approach the problem of Coulomb 
correlations in narrow-band compounds is to divide it in two part: 

(i) mapping of conventional electronic structure calculations onto the multi-orbital Hubbard model, and derivation of 
the parameters of this model from the first principles, for example starting from the simplest electronic structure in 
LDA; 

(ii) solution of this multi-orbital Hubbard model using modern many-body methods. 9 
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In this paper we will discuss the first part of this project and show how results of conventional LDA calculations 
for the t2 g bands can be mapped onto the multi-orbital Hubbard model: 

RR' a/3 R ap-yS 

where c Rct (cr q ) creates (annihilates) an electron in the Wannier orbital Hrq, of the site R, and a is a joint index, 
incorporating all remaining (spin and orbital) degrees of freedom. The matrix ||/irr/|| parameterizes the kinetic energy 
of electrons. The matrix elements fe^, have the following meaning: the site-diagonal part (R=R') describes the 
local level-splitting, caused by the crystal field and (or) the spin-orbit interaction, while the off-diagonal part (R^R') 
stands for the transfer integrals (or the transfer interactions). U a f3js= f dr J dr'W^(r)Wp(r)v SCI (r—r')WMr')Ws(r') 
are the matrix elements of screened Coulomb interaction i> scr (r— r'), which are supposed to be diagonal with respect 
to the site indices. In principle, the off-diagonal elements can be also included into the model. However, we do not 
consider them in the present work. In Sec. IVIIII we will discuss several open questions related with the definition of 
the intersite Coulomb interactions in LDA. 

The first part of this paper will be devoted to derivation of the parameters of the kinetic energy. Then, we will 
explain how to construct the Wannier functions (WFs), which generate these parameters after applying to the KS 
Hamiltonian in the real space. The next part will be devoted to calculations of screened Coulomb interactions, using 
the WF formalism. 



II. THE LMTO METHOD 

In this section we will briefly review the main ideas of linear-muffin-tin-orbital (LMTO) method, as they will be 
widely used in the subsequent sections for the construction of transfer interactions and the Wannier orbitals. The 
method is designed for the solution of KS equations in LDA. For the details and recent developments, the reader is 
referred to the activity of the O. K. Andersen group at the Max-Plank Institute in Stuttgart^ 

Majority of modern electronic structure methods use some basis. The basis functions of the LMTO method, 
(the so-called muffin-tin orbitals - MTOs) have many similarities with orthogonalized atomic orbitals. As we will 
see below, the LMTO method is very convenient for constructing the WFs, and for certain applications, the basis 
function of the LMTO method, from the very beginning, can be chosen as a Wannier function. 

The conventional LMTO approach employs the atomic-spheres-approximation (ASA) , which assumes that the whole 
space of the crystal can be filled by overlapping atomic spheres (Fig.©, so that the overlap between the spheres as 
well as the empty spaces, which are not encircled by any spheres, can be neglected. The MTOs are constructed from 
solutions of KS equations inside atomic spheres (the partial waves), calculated at some energies E v l (typically, the 
center of gravity of the occupied band or of the entire band) , 4>rl , and their energy derivatives <Prl ■ In each atomic 
sphere, the KS potential is spherically averaged. Therefore, the solutions are proportional to the angular harmonics, 
which are specified by the indices L=(£,m) (correspondingly, orbital and azimuthal quantum numbers) . At the 
atomic sphere boundaries, {</>rl} and {</>rl} match continuously and differenciably onto certain envelop functions. 
The latter are typically constructed from irregular solutions of Laplace equation, which rapidly decay in the real space 
(Fig.©. 

It is easy to verify that the functions 0rl and <Prl obey the following "LMTO algebra" : 

(^RL|*Ri) = l, (4) 

(^rlI^rl) = (^rlI^rl) =0, (5) 

(0Rt|^RL) = PRL, (6) 

(H KS - E vL ) |<£ RL ) = 0, (7) 

and 



{H KS - E uL ) |^ Ri ) = |</>R L ), 



FIG. 2: Construction of basis functions of the LMTO method. It is assumed that the whole space of the crystal is filled by the 
atomic spheres. The left panel shows a concrete example of such filling in the VO2 plane of SrV03. The central part of the basis 
function (the "head") is constructed as the solution of Kohn-Sham equation inside atomic sphere R (0rx), which is calculated 
in the center of gravity of the band and whose angular character is denoted by L. The "tail" of the basis function at the 
neighboring site R/ is constructed from the energy-derivatives of 4>WL> (denoted as 4>-r.'l')- At the atomic spheres boundaries 
(shown by Stm and So for the central transition-metal site and neighboring oxygen sites, respectively), 4>rl and <f>u.'L' match 
continuously and differenciably onto certain envelop function (typically, an irregular solution of the Laplace equation). 

where i?KS=^Ks( r ) is the Kohn-Sham Hamiltonian in the real space. Then, one possible choice of the muffin-tin 
orbitals is 

\ X ) = \4>) + \4>){H-E U ). (9) 
In this paper we use the shorthanded notations by Andersen et ai. For instance, Eq. @ should be read as follows: 

\Xrl) = \<hu.) + E I^R'L'KWr'r - Sb/b.Sl'lEul)- 

R'L' 

The first and second terms in the right-hand side of this equation are sometimes called, correspondingly, the "head" 
and the "tail" of MTO. The angular character of MTO at the central site R is specified by that of partial wave </>rl- 
The matrix elements of Hks in the MTOs basis © can be immediately derived by using the properties Q-® °f 
0rl and <p RL : 

(x\Hks\x) =H+(H- E v )E v p{H ~ E u ), 

where E v and p are the diagonal matrices constructed from {E^l} and {prl}, and (xI-^ksIx) is the shorthanded 
notation for the matrix ||(xrl|#ks|xr'L')||- The corresponding overlap matrix, ||(xrz,|xr'L')II> is 

(x\x) =1 + (H- E v )p{H - E v ). 

Since the second term in the right-hand side of (xlx) is typically small, the basis functions © are said to form a nearly 
orthogonal representation of the LMTO method, and H is the LMTO Hamiltonian in the second order of (lhi—E u ). 
The basis © can be orthonormalized numerically, by applying the transformation 

lx)-lx) = lx)<xlxr 1/2 - 

The corresponding LMTO Hamiltonian, 7?=(x|i?Ks|x): which is formally valid in all orders of (H—E^), is given by 

^ = (xixr 1/2 <xi#Ksix>(xixr 1/2 - 

This Hamiltonian will be used as the starting point in the next section, for the definition of transfer interactions 
between certain Wannier orbitals. We start with the formal description of the downfolding method. The construction 
of the Wannier basis functions, underlying this approach, will be considered in Sec. IIVI where we will use again all 
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merits of the LMTO method and show that proper WFs can be constructed by retaining the "heads" of MTOs and 
attaching to them different "tails" . 

For periodic crystals, it is convenient to work in the reciprocal (k) space. Therefore, if it is not specified otherwise, 
we assume that the LMTO Hamiltonian is already constructed in the reciprocal space, after the Fourier transformation 
of MTOs: 

R. 

where N is the number of sites. 

For all considered compounds, results of our ASA-LMTO calculations are in a good agreement the ones obtained 
using more accurate full-potential methods. In our definition of the crystal-field splitting we go beyond the conventional 
ASA and take into account nonsphericity of the electron-ion interactions (see Sec.lvl). 



III. DOWNFOLDING METHOD FOR THE KINETIC-ENERGY PART 



Parameters of the kinetic energy are obtained using the downfolding method, starting from the electronic structure 
in LDA. In order to describe properly the electronic structure of the TM oxides in the valent part of the spectrum 
using the LMTO method, it is typically required several tens or even hundreds basis functions (including the ones 
associated with empty spheres, which are added in order to improve the atomic spheres approximation for loosely 
packed atomic structures). Several examples of such bases will given in Secs. lVI AIVIDI 

What we want to do next is to describe some part of this electronic structure by certain tight-binding (TB) 
Hamiltonian h, which, contrary H, is formulated in the basis of a very limited number of orthogonal atomic-like 
orbitals. For example, in order to reproduce the ti g bands located near the Fermi level, one would like to use only 
three t ^g or bitals centered at each TM site. These orbitals have a meaning of Wannier orbitals, which will be considered 
in Sec. fcl 

We start with the identity by noting that any eigenstate of the LMTO Hamiltonian H can be presented as the sum 
W—\ipt) + \^r), where \ipt) is expanded over the LMTO basis function of the ^g-type, {|xt)} (here, the character of 
the basis function is specified by its "head"), and \ip r ) is expanded over the rest of the basis functions {|Xr)}- Then, 
the matrix equations for LMTO eigenstates can be rearranged identically as 

(H tt - Lu)\ip t ) + H tr \i> r ) = 0, (10) 
H rt \^t) + {H rr ~ uo)\i, r ) = 0. (11) 

By eliminating \ip r ) from Eq. 1)11(1 one obtains the effective w-dependent Hamiltonian in the basis of £2g-states 

fl3r(w) = H u - H tr (H rr - u))- l H rt 

and the "overlap" matrix 

§(w) = 1 + H tr (H rr - u)- 2 H r \ 

satisfying the condition (tp t \S\ipt} = ^- Then, the required TB Hamiltonian, h, is obtained after the orthonormalization 
of the vectors |^> t )— ►|i/>t)= l § 1 ' 2 |^>t) and fixing the energy u> in the center of gravity of the t 2g band (cjq): 

h = S-^MH^MS- 1 / 2 ^). (12) 

Typically, the downfolding is performed in the reciprocal space, and the parameter loq may also depend on k. The 
Hamiltonian hk can be Fourier transformed back to the real space: 

W = £ eik(R ~ R,) ^. 

k 

The site-diagonal part of /irr/ shall describe the crystal-field (CF) splitting caused by the lattice distortion and 
associated with the transfer interactions between t- and r-orbitals, which are eliminated in the downfolding method, 
while the off-diagonal elements have a meaning of transfer interactions. The first application of this approach has 
been considered in Ref. In Sec. I VII we will illustrate abilities of this method for several types of ti g compounds. 

The crystal-field splitting may have another origin, which is related with nonsphericity of the electron-ion 
interactions^ This contribution will be considered in Sec. IVl 
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IV. WANNIER FUNCTIONS 

In the previous section we have shown that there is a TB Hamiltonian, h, formulated in some basis of Wannier 
orbitals {W}. It allows to generate the electronic structure of isolated bands, which is practically identical to the 
electronic structure obtained after the diagonalization of the total LMTO Hamiltonian H . 

In this section we will solve an inverse problem and construct the basis {W}, which after applying to the original 
KS Hamiltonian, generates the matrix h: 

h= (W\H KS \W). (13) 

In a close analogy with the LMTO method, we will first introduce the orbitals {W}, which are related with {W} 
by the orthonormalization transformation 

\W) = \W)(W\W)~ 1/2 , (14) 

and search \ W) in the form: 

N r 

\W) = \W t ) + J2 T r\Wr), (15) 
r=l 

where W t is constructed entirely from the <2g-type solutions of KS equations inside atomic spheres and their energy- 
derivatives {4>t, </>*}: i.e. both 4> t and <pt belong to the TM sites, and 't' stands for the 3d-t2 g partial waves. Each 
W r is constructed from the rest of the partial waves {</>,-, </> r }. Then, Wt and W r can be found from the following 
conditions: 

1. We request only the i-part of \W) to contribute to the matrix elements of the KS Hamiltonian, and search it in 
the form of MTO: 

\W t ) = \^) + \^)(i,-E vt ). (16) 

In this definition, \Wt) is a function of (yet unknown) matrix (), which will be found later. The matrix elements of 
Hks in the basis of these Wannier orbitals are given by: 



(W\H KS \W) = & + (&- E ut )E ut p t (i) - Eut). (17) 

2. The r-parts of the WF, {Wr}, do not contribute to the matrix elements l|17|) . They are introduced only in order 
to make the WFs ff 1 51) orthogonal to the rest of the eigenstates of the Hamiltonian H . Therefore, we search \ W r ) in 
the form: 

\W r ) = \4> r ) +a r \<j) r ), 

where 



1+y/l- AEl rPr 

is obtained from the condition (W r |-ffKs| W r )=0. 

3. The coefficients r r are found from the orthogonality condition of \W) to N r eigenstates {IV'i)} of the original 
LMTO Hamiltonian H: 

N r 

Y J (^\Wr)T r =-mW t ), i=l, iV P . (18) 

r=l 

This allows to include the r-components of the WFs in a systematic way. For example, by taking into consideration 
the 2p-partial waves inside oxygen spheres (iV r =9 wavefunctions for cubic perovskites), the WFs can be orthogonalized 
to N r =9 0(2p) bands, etc. For a given f), the problem is reduced to the solution of the system of linear equations 

CHJ. 

Since \W t ) contributes to Eq. 118|) . the coefficients {r r } will also depend on f). Therefore, the total WF is an 
implicit function of the matrix rj: |VF)=|VF([j)). 

4. The last step is the orthonormalization Ijl4|l . which after substitution into Eq. Ijl3(l yields the following equation 
for the matrix f): 

t) = (wmwm 1 ' 2 h (wmwm 1 ' 2 E vt )E vt p t § - E vt ). 

This equation is solved iteratively with respect to f). 
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A. Spacial extension of Wannier functions 

The choice of the WFs as well as their extension in the real space is not uniquely defined. For many practical 
applications one would like to have "maximally localized" orbitals, 13 although in the context of the WFs, the term 
"maximally localized" itself bears certain arbitrariness and is merely a mathematical construction, because depending 
on the considered physical property one can introduce different criteria of the "maximal localization" . 

Although we do not explicitly employ here any procedure which would pick up the "most localized" representation 
for the WFs, our method well suits this general strategy and the obtained WFs are expected to be well localized 
around the central TM sites. 

There are several ways of controlling the spacial extension of the WFs in the LMTO method. 

1. By using different envelop functions one can, in principle, change the spacial extension of MTOs (Fig. EJ, which 
controls the decay of the original LMTO Hamiltonian in the real space. For example, instead of irregular solutions of 
Laplace equation, one can use Hankel functions of the complex argument. However, any choice should satisfy certain 
criteria of the completeness of the basis set. From this point of view, the use of the Hankel functions is not well 
justified as it typically deteriorates the accuracy of LMTO calculations. Therefore, in the present work we leave this 
problem as it is and fix the LMTO basis set. 

2. Once the LMTO basis is fixed, the relative weight of the TM d-states and other atomic states which contribute 
to the t2 g band cannot be changed (see Fig. For example, the contribution of the oxygen 2p-states cannot be 
replaced by muffin-tin orbitals centered at the TM sites and vise versa. The same proportion of atomic orbitals should 
be preserved in the WFs, constructed for this ti g band. Then, the only parameter which can be controlled is how 
many WFs, centered at different sites of the lattice, contribute to the density of d-states at the given TM site. Then, 
the definition "localized orbital" mean that it is mainly centered around given TM site. Conversely, the "delocalized 
orbital" may have a long tail spreading over other TM sites. Then, it is easy to see that our procedure corresponds to 
the former choice. Indeed, in the first order of (f)— E„ t ) and neglecting for a while the nonorthogonality to the rest of 
the electronic states, the norm of the WF can be obtained from Eq. (|16() as (WlW) — ^^^ , meaning that the WF is 
fully localized at the central TM site. Then, it holds t)=h, which is valid in the second order of (fj— E ut )w& Therefore, 
the leading correction to the above approximation, which define the actual weight of the WF at the neighboring TM 
sites is controlled by the parameters of the kinetic energy h, and is of the order of hpth. As we will see below, the 
latter is small. The conclusion is rather generic and well anticipated for the strongly-correlated systems for which the 
kinetic-energy terms is generally small. 

3. The angular character of the WF at the central TM site should be consistent with the one extracted from the local 
density of states in the region of ti g bands (in the other words, the local density of states at the TM sites should be 
well represented by atomic orbitals {|x*)} used in the downfolding method). Therefore, we choose {|x*)} as the set of 
atomic orbitals which mainly contribute to the local density of states in the region of t% g bands. For these purposes, 
at each TM site we calculate the density matrix in the basis of five d orbitals {|xd)} : 

^=J2<Xd\^i\Xd), (19) 

and sum up the contributions of all ti g bands (here, i is an joint index, which incorporates the band index and the 
coordinates of the k-point in the first Brillouin zone). Then, we diagonalize A/", and assign three most populated 
orbitals, obtained after the diagonalization to {|xt)}- 

V. CRYSTAL-FIELD SPLITTING CAUSED BY NONSPHERICITY OF ELECTRON-ION 

INTERACTIONS 

The contribution of Coulomb interactions to the crystal-field splitting is a tricky issue. Despite an apparent 
simplicity of the problem, one should clearly distinguish different contribution and not to include them twice, in the 
kinetic and Coulomb parts of the model Hamiltonian © . The use of full-potential techniques does not automatically 
guarantee the right answer. However, the atomic-spheres-approximation, which typically supplements the LMTO 
method, will also require additional corrections for the crystal-field splitting. In this section we would like to make 
two comments on this problem. 

1. The nonsphcricity of on-site Coulomb interactions is already included in the second part of the model Hamiltonian 
©. The problem will be discussed in details in Sec. IVIll Therefore, in order to avoid the double counting, the 
corresponding contribution to the kinetic-energy part should be subtracted. From this point of view the use of the 
spherically averaged KS potential in ASA is well justified. 
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The same is true for the intersite Coulomb interactions, if they are explicitly included to the model Hamiltonian 
©■ 

2. All remaining interactions should generally contribute to the crystal-field splitting. In ASA, the proper correction 
at the site R can be found by considering the matrix elements of the Coulomb potential produced by all other atomic 
spheres (or ions) at the site R, 

where Wr=Wr,(i") is the WF centered at the site R, and Z^, is the total charge associated with the sphere R': 
namely, the nuclear charge minus the electronic charge encircled by the atomic sphere. The nonsphcrical part of this 
integral can be easily calculated in the real space by using the multipole expansion for |R + r — R'| _1 . 

In all forthcoming discussions, unless it is specified otherwise, the matrix elements of the crystal-field splitting will 
incorporate the correction A/irr associated with nonsphericity of the electron-ion interactions. 

VI. APPLICATIONS TO TRANSITION-METAL OXIDES 



A. Cubic Perovskites: SrVOa 



SrVC>3 is a rare example of perovskite compounds, which crystallizes in the ideal cubic structure. It attracted a 
considerable attention in the connection with the bandwidth control of the metal- insulator transition. 5 

For the cubic compounds, the separation of the LMTO basis functions into {|xt)} an d {IxV)} used in the downfolding 
method is rather straightforward: three ti g orbitals centered at each V site form the subspace of {|x*)} orbitals, and 
the rest of the basis functions are associated with {|xV)}- 

Parameters of LMTO calculations for SrV03 are given in Tabled The corresponding electronic structure is shown 



TABLE I: Atomic positions (in units of cubic lattice parameter a=3.842A), atomic radii (in A) and basis functions included 
in LMTO calculations for cubic SrVOs. 



type of atom 


position 


atomic radius 


LMTO basis 


number of atoms 


Sr 


(0.5,0.5,0.5) 


1.919 


5s5p4d4/ 


1 


V 


(0,0,0) 


1.470 


4s4p3d 


1 





(0.5,0,0) 


1.032 


2s2p 


3 



in Fig. [3J The downfolding procedure is nearly perfect and well reproduces the behavior of three ti g bands. As 



SrVO, 



V(xy) 



O(x) 




V(xy) 



V(xy) 




R M 



FIG. 3: Left panel: Crystal structure of cubic perovskites and atomic wavefunctions mediating transfer interactions between 
V(t2 9 ) orbitals. The standard V (xy)-0(y)-\ (xy) and V (xy)-0(x)-V (xy) interactions operate in the x- and y-directions, re- 
spectively. The V(xy)-Sr(3z 2 -r 2 )-V(xy) interaction operate in the "forbidden" z-direction. Right panel: LDA energy bands 
for SrV03 obtained in LMTO calculations and after the tight-binding (TB) parametrization using the downfolding method. 
Notations of the high-symmetry points of the Brillouin zone are taken from Ref. 

expected for cubic perovskite compounds^ the transfer interactions between different tig orbitals are small. The 
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dispersion of tig bands is well described in terms of three interaction parameters t\ , t[ , and ti : 

e X j,(k) — 2ii(cos ak x + cosak y ) + 2t\ cosak z + At 2 cosak x cosak y 

(a being the cubic lattice parameter; similar expressions for the yz and zx bands are obtained by cyclic permutations 
of the indices x, y, and z). The parameters t\, t'i, and t 2 , obtained after the Fourier transformation, are listed in 
Table UTI As expected, the nearest-neighbor (NN) dc?7r-interaction t\ mediated by the oxygen 2p-states is the strongest. 



TABLE II: Parameters of transfer interactions for SrVOa (in eV). 



tl 


t[ 




-0.209 


-0.023 


-0.084 



For the xy-orbitals, it operates in the x and y directions. However, there is also an appreciable dd(5-interaction t[ 
operating in the "forbidden" direction (for example, the direction z in the case of xy orbitals). These interactions are 
mediated by the Sr(4d) states and strongly depend on the proximity of the latter to the Fermi level. Therefore, it 
is not quite right to say that the transfer interactions between t 2g orbitals are strictly two-dimensional in the cubic 
latticed Since the La(5ei) states are located even lower in energy than the Sr(4d) ones, the interaction t[ is expected 
to be even stronger in LaTiC>3. However, in the case of LaTiC>3 we have an additional complication associated with 
the orthorhombic distortion. As we will see below, it changes the conventional form of transfer interactions expected 
for the simplified cubic perovskite structure dramatically. 

The corresponding WF is shown in Fig. 21 In this case, the i-part of the WF was constructed from the 3d-t 2g 





FIG. 4: Contour plot of the xy-Wannier orbital of SrV03, in the (001) (left) and (110) (right) planes. The solid and dashed 
line correspond to the positive and negative values of the Wannier function. Atomic positions are shown by symbols. Around 
each site, the Wannier function increases/decreses with the step 0.04 from the values indicated on the graph. 

partial waves inside V spheres. The partial waves of the Sr(4d5s), V(3d-e g ), and 0(2p) types were included into the 
r-part, in order to enforce the orthogonality of the WF to the bands of the aforementioned type. 

Since the tig band is an antibonding combination of the atomic V(3d-t 2g ) and 0(2p) orbitals, the WF has nodes 
located between V and O sites. Fig. [S] illustrates the spacial extension of the WF. It shows the electronic charge 
accumulated around the central V site after adding every new sphere of the neighboring atomic sites. Since the WF 
is normalized, the total charge should be equal to one. In the case of SrVC>3 , 77% of the this charge belongs to the 
central V site, 16% is distributed over four neighboring oxygen sites, about 5% belongs to the next eight Sr sites, and 
1% - to the eight oxygen sites located in the fourth coordination sphere. Other contributions are small. 

It is also instructive to calculate the expectation value of the square of the position operator, (r 2 ) = (TU|r 2 |TU), 
which characterizes the spread of the WF in the method of Marzari and Vanderbilt^ They proposed to define the 
"maximally localized" Wannier orbitals as the ones which minimize (r 2 ). Using the WF shown in Fig. 01 we obtain 
(r 2 )=2.37 A 2 . Unfortunately at present all applications of the method by Marzari and Vanderbilt to the TM oxides 
are limited by MnO^i Therefore, we can make only indirect comparison between two different compounds. The values 
of (r 2 ) reported in Ref . ItH for individual WFs centered at the Mn and O sites were of the order of 0.6-0.8 A 2 , that is 
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FIG. 5: Spacial extension of £2 9 -Wannier function in SrVOs: the electronic charge accumulated around the central V site after 
adding every new sphere of neighboring atomic sites. 



considerably smaller than 2.37 A 2 obtained in our work for SrVC>3. However, such a difference is not surprising. 

1. Our scheme of constructing the WFs is not based on the minimization of (r 2 ). Therefore, our values of (r 2 ) should 
be generally larger. 

2. More importantly, the spacial extension of the WFs depends on the dimensionality of the Hilbert space, which is 
used in the construction of the Hubbard model. For example, we will show in Sec. IVIDl that by treating explicitly 
the e G g states in V2O3 one can easily find more compact representation for the Wannier orbitals. This characteristic 
can be further improved by including the 0(2p) states explicitly into the Wannier basisji^ like in Ref. [jj for MnO. 
However, there is a very high price to pay for this extra localization. This is the dimensionality of the Hilbert space, 
which becomes crucial in the many-body methods for the numerical solution of the model Hamiltonian J3j) . 

Thus, we believe that our WFs constructed for isolated ti g band are indeed well localized. 

For cubic perovskites, there are several ways of extracting parameters of transfer interactions from the first-principles 
electronic structure calculations. For example, one can simply fit the LDA band structure in terms of a small number 
of Slater-Koster interactions^ However, the situation becomes increasingly complicated in materials with the lower 
crystal symmetry, like the orthorhombically distorted perovskite oxides, corundum, or pyrochlore compounds. First, 
the number of possible Slater-Koster interactions increases dramatically. Second, the form of these interactions 
becomes very complicated and differs substantially from the cubic perovskite compounds (one example is the mixing 
of tig and e g orbitals by the orthorhombic distortion, which does not occur in cubic perovskites). Therefore, it seems 
that for complex compounds the only way to proceed is to use the downfolding method. In the next sections we will 
consider several examples along this line. 



B. Orthorhombically Distorted Perovskites: YTiOs 



YTiC>3 is a ferromagnetic insulator. The resent interest to this compound has been spurred by the behavior of 
orbital polarization, which is closely related with the origin of the ferromagnetic ground state. YTiC>3 is typically 
considered in combination with LaTiC>3, which is an antiferromagnetic insulator. The magnetic behavior of these 
two, formally isoelectronic materials, is not fully understood^ 

Contrary to SrVC>3, both YTiC>3 and LaTiC>3 crystallize in the strongly distorted orthorhombic structure (shown 
in Fig. [|>] for YTiC>3, the space group No. 62 in the International Tables; the Schonflies notation is D^). In this 
section we will illustrate abilities of the downfolding method for distorted perovskite compounds, using YTi03 as an 
example. Parameters of LMTO calculations for YTiC>3 are given in Tabic UTTl 

A new problem we have to address here is how to separate the basis functions of the LMTO method onto the 
and {|xr)} orbitals. Note that although the t2 g band is well separated from the rest of the electronic structure also 
in the case of YTi03, the atomic ti g and e g orbitals are strongly mixed by the crystal- field effects and the transfer 
interactions in the distorted perovskite structure. Therefore, the conventional separation into atomic ti g orbitals and 
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FIG. 6: Left panel: Fragment of crystal structure of YTiOa (a "distorted cube"). The inset shows the directions of orthorhombic 
axes. Right panel: LDA energy bands for YTiOa obtained in LMTO calculations and after tight-binding (TB) parametrization 
using the downfolding method. Notations of the high-symmetry points of the Brillouin zone are taken from Ref. IlH 



TABLE III: Atomic positions (in units of orthorhombic lattice parameters a=5.316, 6=5.679, and c=7.61lA), atomic radii (in 
A) and basis functions included in LMTO calculations for orthorhombically distorted YTiC>3. The symbol 'Em' stands for the 



empty spheres. 


type of atom 


position 


atomic radius 


LMTO basis 


number of atoms 


Y 


(0.521,0.073,0.250) 


1.758 


5s5p4rf4/ 


4 


Ti 


(0,0,0) 


1.507 


4s4p3d 


4 





(0.379,0.458,0.250) 


1.052 


2s2p 


4 





(0.309,0.191,0.558) 


1.043 


2s2p 


8 


Em 


(0.192,0.268,0.116) 


0.595 


Is 


8 


Em 


(0.438,0.329,0.750) 


0.757 


ls2p 


4 


Em 


(0.149,0.368,0.027) 


0.556 


Is 


8 


Em 


(0.106,0.208,0.199) 


0.479 


Is 


8 



the rest of the basis functions does not apply here, and in order to generate {|x*)} we use eigenvectors of the density 
matrix (see Sec. II V AT) . 

For the site 1, shown in Fig.® these three u t2 g orbitals" have the following form (in the basis of \xy), \yz), \z 2 ), 
\zx), and \x 2 —y 2 ) orbitals, in the orthorhombic coordinate frame): 

|xi) = ( 0.13,-0.60,-0.24, 0.34,-0.67), 

|X2) = ( 0.17, 0.50,-0.35, 0.77, 0.11), (20) 
| X3 ) = (-0.43,-0.54,-0.29, 0.22, 0.62). 

At the sites 2,3, and 4 similar orbitals can be generated from the ones at the site 1 using the symmetry properties of 
the group and applying the 180° rotations around the orthorhombic axes a, c, and b, respectively. These orbitals 
define the local basis (or the local coordinate frame) around each Ti site. 

Then, the rest of the basis functions {x r } can be eliminated using the downfolding method. The corresponding 
electronic structure for the tig bands is shown in Fig. which reveals an excellent agreement between results of 
LMTO calculations and their tight-binding parametrization using the down-folding method. 

Parameters obtained after the transformation to the real space are listed in Table llVl in the local coordinate frame. 
We note a substantial crystal-field splitting associated with the orthorhombic distortion. After the diagonalization 
of the site-diagonal part of the TB Hamiltonian h, we obtain the following ("one-down, two- up") splitting of ti g 
levels: —0.076, 0.032, and 0.046 eV. Some implications of the crystal-field splitting to the orbital polarization and 
the magnetic ground state of YTiC>3 and LaTiOa have been discussed in Refs. [Tl| and IlIi 1 ^ The form of transfer 
interactions becomes extremely complicated, and differs dramatically from many naive expectations based on the 
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TABLE IV: Parameters of the kinetic energy for YTiC>3 (in eV), in the local coordinate frame, the atomic positions are shown 
in Fig. E] The basis functions at the site 1 are given by Eqs. H2()fl. The basis functions at the sites 2 and 3 are obtained by 
the 180° rotations of the site 1 around the orthorhombic axes a and c, respectively. The matrix hn describes the crystal-field 
splitting at the site 1 (nonsphericity of electron-ion interactions has been added in /in) The matrices hi2 and hi3 stand for the 
transfer integrals in the bonds 1-2 and 1-3, respectively. 



i-j 




fin 






-0.060 


-0.039 


-0.014 


1-1 


-0.039 


0.028 


-0.001 




-0.014 


-0.001 


0.032 




-0.080 


-0.036 


0.072 


1-2 


0.162 


0.013 


0.073 




0.063 


0.031 


-0.038 




0.000 


-0.080 


0.055 


1-3 


-0.080 


0.082 


-0.007 




0.055 


-0.007 


0.073 



analogy with the cubic perovskites. Generally, the transfer interactions are three-dimensional and operate between 
different ti g orbitals. 

The WFs are shown in Fig. and their spacial extension in the real space is illustrated in Fig. [5] In these 
calculations, the WFs have been orthogonalized to the 0(2p), Ti(3d-e ff ), and Y(4d) bands. The orbitals appear 




'---0.02 



FIG. 7: Contour plot of three t2 9 -Wannier orbitals in the (001) plane of YTiOs. The orbitals are centered at the site 1 shown 
in Fig. [5] Projections of the vanadium and oxygen sites onto the plane are shown by circles and crosses, respectively. Other 
notations are the same as in Fig. 2] 

to be more localized than in SrVOa: 82-87% of the total charge belongs the central Ti site, and only 6 to 10% is 
distributed over neighboring oxygen sites. This is because of the large orthorhombic distortion, which suppresses all 
interatomic interactions mediated by the oxygen states. Another reason is the larger energy distance between 0(2p) 
and t2 g bands in YTi03 (3.2 eV against 0.3 eV in SrV03 - see Fig.QJ, which explains smaller weight of the atomic 
oxygen states in the t^g band and the WFs of YTiOs. 

Another interesting feature is that the degree of localization is pretty different for three orbitals. For example, we 
obtain (r 2 )= 2.28, 1.90, and 2.05 A 2 , correspondingly for Wi, W2, and W3 shown in Fig. One can paraphrase it 
in a different way: the degree of hybridization can be different for different ti g orbitals, unless they are related with 
each other by symmetry operations. 

C. Pyrochlores: Y 2 Mo 2 7 

The pyrochlore compounds exhibit a variety of interesting properties. Many of them are not fully understood. 
Y2M02O7 is a canonical example of geometrically frustrated systems. In this compound, the magnetic atoms form 
the networks of corner-sharing tetrahedra. Then, the antiferromagnetic coupling between NN Mo spins leads to the 



13 




FIG. 8: Spacial extension of t2 9 -Wannier functions in YTiOa. 



frustration. The origin of this antifcrromagnctic coupling can be understood on the basis of semi-empirical Hartree- 
Fock calculations^ A remaining question, which is not fully understood, is the origin of the spin-glass state realized 
in Y2M02O7 below 20 K. ? A For comparison, Nd2Mo2C>7 is a ferromagnet, revealing a large anomalous Hall effecti22i22 
Another interesting group is superconducting /3-pyrochlores with the chemical formula ^40s2C>6 (A= K, Rb, and 
Cs)& 

In this section we will derive parameters of the kinetic energy for the ti g band of Y2M02O7. Very similar strategy 
can be applied for other pyrochlores. Parameters of LMTO calculations for Y2M02O7 are listed in Table Ivl 

TABLE V: Atomic positions (in units of cubic lattice parameter a=10.21 A), atomic radii (in A) and basis functions included 
in LMTO calculations for pyrochlore Y2M02O7. 



type of atom position atomic radius LMTO basis number of atoms 

Y (0.5,0.5,0.5) 1.746 5s5p4d 4 

Mo (0,0,0) 1.587 5s5p4d 4 

O (0.375,0.375,0.375) 1.032 2s2p 2 

O (0.338,0.125,0.125) 1.032 2s2p 12 

Em (0.125,0.125,0.125) 0.952 Is2p3d 2 

Em (0.25,0,0) 1.044 Is2p3d 8 



In the pyrochlore lattice, each Mo site is located in the trigonal environment (Fig.|5J. Therefore, the atomic Mo(t2 9 ) 
levels will be split into one-dimensional a\ g and two-dimensional representations. The latter states can mix with 
the Mo(e^) states, which belong to the same representation. Therefore, the basis functions {|xt)} can be constructed 
in the same manner as for YTiC>3, by diagonalizing the site-diagonal part of the density matrix for the ti g bands. For 
the site 1 shown in Fig. this yields the following atomic orbitals (in the basis of \xy), \yz), |z 2 ), \zx), and \x 2 —y 2 ) 
orbitals): 



Xi) 


= ( 


0.58, 


0.58, 


, 


0.58, 


), 




= (- 


-0.06, 


0.18, 


-0.29, 


-0.13, 


0.93), 


£3) 


= ( 


0.18, 


-0.04, 


0.93, 


-0.14, 


0.29). 



In these notations, the first orbital correspond to the a\ g representation, and two other - to the e g representation. 
Similar orbitals at the sites 2, 3, and 4 can be generated from the ones at the site 1 using the symmetry operations 
of the 7 h group (No. 227 in the International Tables): namely, the 180° rotations around the cubic axes x, y, and z, 
respectively. The rest of the basis functions form the subspace {|xV)}- 

The electronic structure obtained after the elimination of the {|Xr)} orbitals is shown in Fig. |5J Again, we note an 
excellent agreement with the results of the original LMTO calculations. The parameters of the kinetic energy in the 
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FIG. 9: Left panel shows the fragment of crystal structure of Y2M02O7: a single Mo-tetrahedron surrounded by the oxygen 
sites. In Y2M02O7, the Y and Mo ions from two sublattices of corner-sharing tetrahedra. The sublattices interpenetrate in a 
way, in which each Mo is located in the center of Y hexagon (not shown here). Right panel: LDA energy bands for Y2M02O7 
obtained in LMTO calculations and after tight-binding (TB) parametrization using the downfolding method. Notations of the 
high-symmetry points of the Brillouin zone are taken from Ref. El 

real space are listed in Table fVTI We note an appreciable (^0.25 eV) crystal-field splitting between the a\ g and e g 

TABLE VI: Parameters of the kinetic energy for Y2M02O7 (in eV), in the local coordinate frame. The atomic positions are 
shown in Fig. El The basis orbitals at the site 1 are given by Eqs. I2H . The basis orbitals at the sites 2, 3, and 4 are obtained 
by the 180° rotations of the site 1 around the cubic axes x, y, and z, respectively. The matrix hn describes the crystal-field 
splitting at the site 1 (nonsphericity of electron-ion interactions has been added in fen). The matrix hi4 stands for transfer 
integrals between atoms 1 and 4. 



i-j 










-0.177 


0.000 


0.000 


1-1 


0.000 


0.088 


0.000 




0.000 


0.000 


0.088 




-0.064 


-0.001 


0.002 


1-1 


-0.001 


-0.186 


-0.070 




0.002 


-0.070 


0.022 



orbitals. 20 The NN transfer interactions in the bonds other than 1-4 can be obtained using the symmetry operations 
of the 7 h group. The transfer interactions beyond the nearest neighbors are considerably smaller. 

The corresponding WFs are shown in Fig. ^3 and their spacial extension is depicted in Fig. ^2 The WFs have 
been orthogonalized to the neighboring Y(4d-e g ) and 0(2p) bands. Since the 4d-wavefunctions are typically more 
extended in comparison with the 3d ones, the WFs are less localized. In the case of Y2M02O7, 75-80% of the total 
charge is located at the central site, and about 20% is distributed over neighboring oxygen sites. The degree of 
localization also depends on the symmetry of WFs. So, the a\ g orbital is well localized within the MoC>6 cluster, 
whereas the e g orbitals have a noticeable weight (~2.5% of the total charge) at the Y and Mo sites belonging to the 
next coordination sphere. 

D. Corundum- type V2O3 

V2O3 is regarded as the canonical Mott-Hubbard system, where the Coulomb interaction between conduction 
electrons leads to a breakdown of the conventional one-electron band theory^ It was and continues to be the subject 
of vast research activity, which has been summarized in many review articles (for instance, Ref. 0). 

V2O3 crystallizes in the corundum structure with two formula units per rhombohedral cell (the space group is D® d , 
No. 167 in the International Tables). The local environment of the V sites is trigonal (Fig. ll2|l . in which the tig levels 
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FIG. 10: Contour plot of three t2 9 -Wannier orbitals in the (001) plane of Y2M02O7. The orbitals correspond to the site 1 in 
Fig. [5] The orbital 1 corresponds to the a\ g representation, the orbitals 2 and 3 correspond to the e g representation. Other 
notations are the same is in Fig. [I] 
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FIG. 11: Spacial extension of f2 9 -Wannier functions in Y2M02O7. 



are split into one-dimensional a\ g and two dimensional e" g representations. 
Parameters of LMTO calculations for V2O3 are given in Table IVTfl 

TABLE VII: Atomic positions (in units of rhombohedral lattice parameters a=2.859 and c=4.667 A), atomic radii (in A) and 
basis functions included in LMTO calculations for V2O3. 



type of atom 


position 


atomic radius 


LMTO basis 


number of atoms 


V 


(0,0,0.5) 


1.317 


4s4p3d 


4 


O 


(-0.5,0.326,0.25) 


1.009 


2s2p3d 


6 


Em 


(0,0,0) 


1.223 


Is2p3d 


2 


Em 


(0.5,0.267,-0.25) 


0.921 


ls2p 


6 



In the present work, our main interest in V2O3 will be purely academic. As we can see in Fig. ^ V2O3 has two 
well separated bands, which are mainly formed by the V(3d) states. One is the ti g band, which in LDA is crossed by 
the Fermi level. Another one is the e g band, which is located around 3 eV, and composed mainly of the states. 
The latter can mix with the e g ones. Therefore, for V2O3 (and related corundum-type oxides) one can introduce two 
different models. 




FIG. 12: Left panel: Crystal structure of V2O3 with the notation of atomic sites. The sites 3 and 4 are equivalent and can be 
transformed to each other by primitive translations (both atomic positions are required in order to specify the form of transfer 
integrals in Table IVim . Right panel: LDA energy bands for V2O3 obtained in LMTO calculations and after tight-binding 
(TB) parametrization using the downfolding method for the three-orbital (top) and five-orbital (bottom) models. Notations of 
the high-symmetry points of the Brillouin zone are taken from Ref. IhH 



The first one is more general and explicitly treats all V(3d) bands. In the following we will call it as the "five- 
orbital" model, according to the number of basis orbitals {|xt)} P er one V site. The oxygen degrees of freedom will 
be eliminated using the downfolding method. The corresponding TB Hamiltonian will be denoted as h^ 5 \ 

The second one is the minimal model, which can be derived from the previous one by eliminating the e g states. We 
will call it the "three-orbital model" . The corresponding TB Hamiltonian will be denoted as . 

The basic difference between these two models is that the first one treats the e g states explicitly, while in the second 
case the effect of these states is included implicitly, through the renormalization of interaction parameters between 
a\ g and e g orbitals. 

For both models, the local orbitals at each V site were obtained from the diagonalization of the density matrix, 
which sums up the contributions over twelve t2 g bands. For the three-orbital model, such choice of the basis functions 
is very important, as it controls the accuracy of the downfolding method. For the five-orbital model, one can use 
any unitary transformation of the five 3d orbitals. Obviously, the final result will not depend on this transformation. 
However, for a better comparison between two models, we use the same basis in both cases. In principle, the e g and 
e g states will be mixed in the density matrix, as they belong to the same representation. However, we will continue 
to call the lower- and upper- lying e g states as e g and e gl despite the fact that each of them may have an admixture 
of another type. 

Then, the basis orbitals at the sites 1 and 2 (see Fig. I12f) have the following form (in the basis of atomic \xy), \yz), 
\z 2 ), \zx), and \x 2 —y 2 } orbitals): 



|Xi) = ( 0.77,-0.29, , 0.48,-0.30), 
|xa) = ( 0.30, 0.48, , 0.29, 0.77), 
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|X3> = ( ,0 ,1 ,0 ,0 ), (22) 
|X4) = ( 0.49, 0.51, ,-0.65,-0.27), 
|X5) = ( 0.27,-0.65, ,-0.51, 0.49). 

In these notations, |xi) and \x_2) are the e g orbitals, |x3) is the a\ g orbital, and |x4) and |xs) are the e g orbitals. The 
basis orbitals at the sites 3, 4, and 5 are generated by the mirror-reflection y^— y in Eq. 1|22|) . In the three-orbital 
model, first three orbitals constitute the subspace while two remaining orbitals arc included in {|xV)}- In the 

five-orbital model, all five orbitals are included in 

The electronic structure obtained after the downfolding is shown in Fig.^| The three-orbital model well reproduces 
the behavior of twelve ti g bands, while the five-orbital model allows to reproduce both ti g and e g bands. 

The corresponding parameters in the real space are given in Table I Villi In the five-orbital model, one can see a 
noticeable hybridization between t 2g and e g states. Therefore, the elimination of the e g states in the three-orbital 
model should lead to an additional renormalization of the parameters of the crystal-field splitting and the transfer 
interactions. Generally, the matrix elements of the kinetic-energy part in the subspace oit2 g orbitals are not the same 
for two considered models. 



TABLE VIII: Parameters of the kinetic energy (in eV) in the local coordinate frame for three- and five-orbital models of V2O3 
(denoted as and ft' 5 ', respectively). The atomic positions are shown in Fig. 1121 The basis functions at the sites 1 and 2 
are given by Eqs. 1221 . The basis functions at the sites 3, 4, and 5 are obtained 1 by the refection y^— y of the site 1. The 



order of orbitals is e^, 


ai g , and 


e%. The 


orbitals 


are eliminated in the three-orbital model. 










ij 






1.(3) 






£(5) 








-0.053 


0.000 


0.000 


0.141 


0.000 


0.000 


0.008 


0.056 




0.000 


-0.053 


0.000 


0.000 


0.141 


0.000 


-0.056 


0.008 


1-1 


0.000 


0.000 


0.107 


0.000 


0.000 


0.190 


0.000 


0.000 










0.008 - 


-0.056 


0.000 


2.326 


0.000 










0.056 


0.008 


0.000 


0.000 


2.326 




0.024 


0.089 


-0.103 


0.009 


0.080 


-0.104 


0.026 


-0.043 




0.089 


-0.150 


0.207 


0.080 - 


-0.103 


0.202 


-0.085 


0.102 


1-2 


-0.103 


0.207 


-0.025 


-0.104 


0.202 


-0.032 


0.114 


-0.068 










0.026 - 


-0.085 


0.114 


-0.067 


-0.030 










-0.043 


0.102 


-0.068 


-0.030 


-0.042 




-0.037 


-0.027 


0.000 


-0.032 - 


-0.054 


0.000 


0.064 


0.163 




-0.027 


0.037 


0.000 


-0.054 


0.032 


0.000 


0.163 


-0.064 


1-3 


0.000 


0.000 


-0.342 


0.000 


0.000 


-0.331 


0.000 


0.000 










0.064 


0.163 


0.000 


-0.014 


-0.065 










0.163 - 


-0.064 


0.000 


-0.065 


0.014 




-0.089 


0.024 


-0.102 


-0.112 


0.013 


-0.113 


-0.053 


-0.002 




0.024 


-0.018 


0.024 


0.013 


0.000 


0.048 


-0.356 


-0.029 


1-4 


-0.102 


0.024 


-0.125 


-0.113 


0.048 


-0.124 


-0.092 


-0.004 










-0.053 - 


-0.356 


-0.092 


-0.213 


-0.031 










-0.002 - 


-0.029 


-0.004 


-0.031 


0.000 




0.053 


-0.053 


0.019 


0.046 - 


-0.060 


0.027 


-0.076 


-0.175 




0.004 


-0.039 


0.034 


-0.002 - 


-0.050 


0.046 


0.069 


0.151 


1-5 


-0.047 


0.113 


-0.054 


-0.071 


0.105 


-0.050 


-0.020 


-0.032 










-0.124 


0.017 


0.055 


-0.039 


-0.078 










-0.259 


0.041 


0.161 


-0.133 


-0.234 



As an example, we shown in Table llXl the crystal-filed splitting between and a\ g levels for the series of corundum- 
type oxides, obtained after the diagonalization of the matrices fif^ and . The splitting turns out to be very different 
in two different models. Since many properties of TM oxides are controlled by this crystal-field splitting^ such a 
model-dependence may be viewed as somewhat unphysical. However, the crystal-field splitting cannot be considered 
independently from other model parameters, such as the Coulomb and transfer interactions, which should be defined 
on the same footing and for the same type of model. For the Coulomb interactions, it is important to follow the 
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concept of WFs, which we will consider in the next section. 

TABLE IX: Parameters of crystal-field splitting between e g and ai g levels obtained for the series of corundum- type oxides in 
the three- (A^p) and five- (Aq f ) orbital models. The positive value of Acp means that the a\ g level lies higher than e g . The 
nonsphericity of electron-ion interactions has been included (the values in parenthesis show the crystal-field splitting originating 
from the transfer interactions alone). 



compound 


Ag> (cV) 


Ag> (eV) 


Ti 2 3 


0.070 (0.254) 


-0.033 (0.159) 


V 2 3 


0.160 (0.233) 


0.051 (0.128) 


Cr 2 3 


0.140 (0.216) 


0.051 (0.129) 



We also note an appreciable contribution coming from nonsphericity of the electron-ion interactions. This contri- 
bution, which is ignored in conventional ASA, acts against the crystal-field splitting originating from the transfer 
interactions and tends to stabilize the ai g level. This may revise certain conclusions obtained in the framework of 
ASA-LMTO method^ 

Corresponding WFs are shown in Fig. ^3 and their spacial extension - in Fig. [21 All functions have been 




FIG. 13: Contour plot Wannier orbitals in the (100) plane of V2O3, obtained in the three-orbital (top) and five-orbital 
(bottom) models. Projections of the vanadium and oxygen sites onto the plane are shown by circles and crosses, respectively. 
The numbering of V atoms ("V 1" and "V 3") is the same as in Fig. 1121 First two orbitals (1 and 2) are of the e g type, third 
orbital (3) is of the a\ g type, and last two orbitals (4 and 5) are of the e g type. Other notations are the same is in Fig. [7| 

orthogonalized to the 0(2p) band. In the three-orbital model, the WFs have been additionally orthogonalized to the 
V(e 9 ) band. Similar to Y2M02O7, the spacial extension of the WFs strongly depends on their symmetry. Generally, 
the ai g and orbitals are more localized, while the e a g orbitals have a considerable weight (more than 20% of the 
total charge) at the neighboring oxygen sites. Furthermore, the spacial extension of the WFs depends on the model 
for which they are constructed. Generally, the five-orbital model allows to construct more compact WFs rather than 
the three-orbital one. For example, in the three-orbital model we have (r 2 )= 1.75 and 2.38 A 2 , correspondingly for 
the a\ g and orbitals. For comparison, the five-orbital yields (r 2 ) = 1.04, 1.04, and 1.41 A 2 for the ai g , e g , and e g 
orbitals, respectively. 
This is not surprising. 

1. The transfer interactions in the three-orbital model are longer-ranged, as they contain additional contributions 
mediated by the e a g orbitals. 

2. For the three-orbital model, the WFs should be additionally orthogonalized to the V(e 9 ) band. At the central 
site, this condition can be easily satisfied by choosing proper atomic t2 g and e g orbitals, which diagonalize the density 
matrix. However, the WF has a tail spreading to the neighboring sites, which should be additionally orthogonalized 
to the V(e 9 ) band by including partial waves of the type into the r-part of the WF. 

Thus, there is certain compromise with the choice of the suitable model for compounds like V2O3, where smaller 
dimensionality of the Hilbert space in the three-orbital model is counterbalanced by necessity to deal with more 
extended WFs. 
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FIG. 14: Spacial extension of Wannier functions in the three- (top) and five- (bottom) orbital models for V2O3. 

VII. THE EFFECTIVE COULOMB INTERACTION 

The calculation of effective Coulomb interactions for the first principles is an extremely complicated problem because 
they are subjected to different mechanisms of screening which should be taken into consideration in the process of 
these calculations. So far, the solution of this problem has not been fully accomplished by any of the research groups, 
despite a vast activity in this directioni 2 ^i 2 Ti 2 ?i 2 ?i?9]? 1 i? 2 ]??i? 4 ]?^i??]?7i??i??i 4 ' i 1 It would be probably fair to say from the 
very beginning that we were not able to solve this problem either, without additional approximations, which will be 
considered in Sec. IVII Al However, we hope to present certain systematics on different points of view, which can be 
found in the literature. We will also summarize several open questions and unresolved problems. 

It is convenient to start with the basic definition of the effective Coulomb interaction U, as it was discussed in 
many details by Herring. 41 According to this definition, (the spherically averaged part of) U is nothing but the energy 
cost for the reaction 2(d n )^d n+1 +d n ~ 1 , i.e. for moving a d-electron between two atoms, located at R and R', and 
initially populated by n R =n R i=n electrons: 

U = E[n R + 1, n w - 1] - E[n R , n R ,}. (23) 

This U may depend on R and R', and using several combinations of R and R' one can extract the values of both 
on-site and intersite interactions. A typical example for SrV03 will be considered in Sec. IVII Bl However, here we 
drop these atomic indices and consider more general aspects of calculations of the effective Coulomb interactions. 
It is implied that the electron is transferred between two Wannier orbitals, and n R and n R > are the populations 
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of these Wannier orbitals. A special precaution should be taken in order to avoid the double counting of the kinetic 
energy term. Indeed, since the kinetic-energy term is included explicitly into the Hubbard model it should not 
contribute to the total energy difference 1)23(1 . This point was emphasized by Gunnarsson and co-workers, in the 
series of publicationsi22*2ii22 They proposed to derive U from constraint-LDA (c-LDA) calculations and suppress 
all matrix elements of hybridization involving the atomic <i-states. Such a procedure can be easily implemented in the 
LMTO method. For the 3d-compounds, this method typically yields C/~5-12 eVj&2£, which is too large (if correct, 
results of this approach would imply that all nature surrounding us would be "strongly correlated"). Therefore, 
although the basic strategy is correct, there is an important piece of physics, which is missing in the method of 
Gunnarsson et al. 

Similar strategy can be pursued in our WF method. Our basic idea is to switch off the kinetic-energy term during 
the construction of the WFs, and to use these functions in calculations of the effective interaction U. Therefore, instead 
of regular WFs {W(h)}, which after applying to the KS Hamiltonian generate the matrix h = (W (h)\Hics\W (h)) , 
we introduce the set of auxiliary Wannier functions {W(c)}, satisfying the condition (W (c)\Hks\W '(c)) = c. In the 
ground-state configuration (nR=riR/=n), c is a constant, which can be dropped^ In the excited state (riR^nR'), 
c is a diagonal matrix with respect to the site indices, c = ||cr <5rr'||, where each matrix element cr may depend 
on occupation numbers {hr}. Since such auxiliary WFs do not interact with each other through the kinetic-energy 
term, they can be used as the basis functions for the effective Coulomb interaction U. 

The auxiliary WFs can be easily constructed using the method proposed in Sec. IIVI after the substitution h=c in 
all equations. Meanwhile, the orthogonality condition to other LDA bands is strictly observed by including proper 
solutions of KS equations inside atomic spheres and their energy derivatives into the r-part of auxiliary WFs. This 
allows to retain the hybridization between TM d- and oxygen p-states, which is an important feature of TM oxides. 
As we will see below, the change of this hybridization, induced by the reaction 2(d n )^d n+1 +d n ~ 1 , represents a very 
important channel of screening, which substantially reduces U and explains many details of its behavior in solids. 
This channel of screening has been overlooked by Gunnarsson et al. 

A similar idea, although formulated in the very different way, has been recently proposed by Aryasetiawan et al~ 
They proposed to extract the parameter U from the GW method, 43 - 44 and suppressed all contributions to the GW 
polarization function associated with the transitions between Hubbard (in our case - t^g) bands, in order to avoid 
the double counting of these effects in the process of solution of the Hubbard model. Clearly, since the polarization 
function in the GW method will vanish without the kinetic-energy term, this procedure appears to be similar to the 
setting h=c for the WFs, which are used as the basis functions for the effective Coulomb interaction U. 

A characteristic example of the auxiliary WFs is shown in Fig. ^] for SrV03. We note only a minor difference 



FIG. 15: An example of auxiliary Wannier function {W'(c)} for SrVC>3 for which all parameters of the kinetic energy are set 
to be zero. Left panel shows the contour plot of the ly-orbital in the (001) plane. Right panel shows the spacial extension of 
this orbital. Other notations are the same as in Figs. 0] and 

between auxiliary WFs and the regular ones shown in Fig. 0] meaning that the main details of the WFs for the t^ g 
bands are predetermined by orthogonality condition to other bands. For example, in the case of auxiliary WFs, 79% 
of the total charge are accumulated at the central V site (instead of 77% for the regular WFs). The values of (r 2 ) 
obtained for the auxiliary and regular WFs are 2.27 and 2.37 A, respectively. Thus, the auxiliary WFs appear to be 
more localized. However, the difference is small. The result is well anticipated for strongly correlated systems, for 
which the kinetic-energy term h is expected to be small. 
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Since the KS Hamiltonian is diagonal in the basis of auxiliary WFs, the latter can be regarded as eigestates of this 
Hamiltonian corresponding to certain boundary condition. The corresponding set of KS eigenvalues will be denoted 
as {er}. Then, the occupation numbers {np.} become well defined and one can use the standard properties of the 
density-functional theory. Namely, by using Janak's theorem for the KS eigenvalues 

SE 

£r 



and Slater's transition-state arguments, Eq. (|23|l can be further rearranged as 

1 1. r 1 1. 

U = SR[n R + 2' n R-' _ 2 £R ^ R - 2'" R ' + 2 

The final expression for the parameter U is obtained by considering the deviations ±1/2 from zjr and nR' as a weak 
perturbation and employing the Taylor expansion. Then, in the first order of ±1/2 one obtains: 

U=P±, (24) 

where the energy derivative is calculated under the following condition: 

^r + Br' = const, (25) 

which guarantees the conservation of the total number of particles. Strictly speaking, the definition (|24|l corresponds 
to the infinitesimal change of the occupation numbers 2(d n )^d n+Sn +d n ~ 5n , which is different from original Herring's 
definition (|2*5|) . However, for practical purposes, these definitions can be regarded as equivalent as they yield very 
similar values for the parameter U ~ 

Then, it is convenient to use the Hellman-Feinman theorem, which allows to relate U with the change of the Hartree 
potential (the change of the exchange-correlation potential in LDA is typically small and can be neglected)^ 

u = m^m. 

Taking into account that Vn(r)=e 2 J dr'p(r')/|r — r'|, and using Eq. (|5J) for the electron density, the above expression 
can be rearranged as 



5n 



r 



Sp(r) f Sni | / f \\2 & i / m2 



l|^(r)| 2 +n i7 ^|^(r)| 2 . (26) 



where 

The last expression points out at the existence of two additive channels of screening. 

(i) The first one comes from the change of occupation numbers. Due to the constraint (|25|l imposed on the occupation 
numbers, this channel involves two Wannier orbitals, centered at different TM sites, and describes the screening of 
on-site Coulomb interactions by intersite interactions. Other states can be affected by this term only through the 
change of wavefunctions in the process of iterative solution of the KS equations. We also would like to note that this 
channel of screening is absent in the GW methods, which may lead to an error for metallic compounds^ 

(ii) The second channel describes the relaxation of the wavefunctions. It affects both the auxiliary WFs and the 
electronic states belonging to the rest of the spectrum. 



A. Approximations and Simplifications 

The usual way in calculating the parameter U is the c-LDA approach, that is to solve iteratively Eqs. and 
(J2J) for a fixed set of occupation numbers {ni}, which does not necessarily follow the Fermi-Dirac distribution for 
the ground s ta,tem£i2L2fL2iL2iL2LMi2& In practical calculations, these occupation numbers are controlled by an external 
potential 5V C xt (r) , playing a role of Lagrange multipliers in the constrained density functional theory. In spite of many 
limitations for the strongly correlated systems, LDA is formulated as the ground-state theory. Therefore, there is a 
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general belief that it should provide a good estimate for the total energy difference given by Eq. Il2.'j[| and all other 
expressions which can be derived from Eq. (|23|l using usual arguments of DFT. 

From the practical point of view, the basic difficulty of combining c-LDA with the WF method is the necessity 
to deal with relaxation of these WFs. This means that the auxiliary WFs should be recalculated on each iteration, 
for every new value of the electron density and the KS potential. Taking into account an arbitrariness with the 
choice of the WFs, this procedure cannot be easily implemented in the standard c-LDA calculations. Instead, we will 
employ a hybrid method, which starts with c-LDA and then takes into account the effects of relaxation of the WF 
in an analytical form, using the well-known expressions for the screened Coulomb interaction in the random-phase 
approximation (RPA). 37 i 38 i 39 i 43 i 44 

In c-LDA calculations, the change of the occupation numbers {Srii} is associated with some change of the total 
potential SV=SV ex t+SVa+SVxc- Then, the change of the electron density in Eq. can be identically expressed in 
terms of the polarization function as 

6p(r) = J dr'P(r,r')5V{r'). (27) 

Using Eq. 12{j|> . one can identify three main contributions to the polarization function associated with the following 
processes (correspondingly P 1 , P 11 , and P 111 ): 

I. the change of the occupation numbers of the auxiliary WFs; 

II. the relaxation of the auxiliary WFs; 

III. the relaxation of the rest of the electronic states. 

Then, each W can be expressed in terms of the basis functions (or partial waves) {</>} and {(/>}, using Eq. JTSJ. 
Therefore, the change of W includes the relaxation of these basis functions as well as the change of hybridization 
of the TM tig states with (mainly) the oxygen states. The latter is given by the change of coefficients {r r } in the 
right-hand side of Eq. (|15|) . The corresponding contributions to the polarization function are denoted as P IIB and 
P IIH ( which stand for the change of basis functions and hybridization, respectively. The same arguments are applied 
to relaxation of the rest of the electronic states. The corresponding polarization function can be divided accordingly 
in P IIIB and P IIIH . We also introduce combined notations: pB—piiB +p uiB and p H =p iiH +p niH 

We would like to point out here that, conceptually, the RPA approach for treating the relaxation effects is similar 
to c-LDA. The main difference is that RPA is based on an analytical expression for the change of the wavefunctions, 
formulated in terms of the perturbation theory expansion, while c-LDA treats the same effects numerically^ There- 
fore, we use a hybrid c-LDA+RPA scheme, which was originally considered in Ref. |4(J. It consists of two steps, 
(i) First, we take into account the screening associated with P 1 , and P B in the framework of conventional c-LDA 
method, and neglect all kinds of hybridization effects involving the TM d-orbitals. An example of such a model elec- 
tronic structure is shown in Fig.^J This part is totally equivalent to the method of Gunnarsson and co-workers 
It allows to calculate the Coulomb repulsion u and the intra-atomic exchange (Hund's rule) coupling j——25 2 E/5m 2 
in the atomic limit (m being the spin magnetization) . By using these u and j one can construct the full 5x5x5x5 
matrix u of Coulomb interactions between atomic d electrons, as it is typically done in the LDA+U methodic This 
matrix will be used as the starting point in RPA calculations. 

c-LDA is supplemented with additional approximations, such as the atomic-spheres approximation. It also disre- 
gards some hybridization effects. However, in Sec. IVIlO we will see that at least for the static Coulomb interaction, 
the RPA results are close to the strong-coupling regime. In such a situation, the precise value of the parameter u, 
which is used as the starting point for these calculations, appears to be less important, and it is sufficient to have an 
"order of magnitude" estimate, which can be obtained from c-LDA. (ii) We turn on the hybridization, and evaluate 
the screening associated with the last portion of the polarization function, P H , in RPA: 



U 



1-uP 



H 1 



u, (28) 



where P H is the 5x5x5x5 matrix P H =\\P^ jS \\, which will be specified below. 

Since total P is an additive function of P 1 , P B , and P H , this procedure can be justified within RPA, where each new 
contribution to the polarization function (P H ) can be included consequently by starting with the Coulomb interaction 
u, which already incorporates the effects of other terms (P 1 and P B )mtL 

The physical meaning of processes associated with the change of the hybridization and their role in the screening 
of local Coulomb interactions is illustrated schematically in Fig. 1171 

Since the creation and the annihilation of an electron in RPA are treated as two independent processes, the screening 
of Coulomb interactions will be generally different from that associated with the true reaction 2(d n )^d n+1 +d n ~ 1 . 




FIG. 16: An example ol model electronic structure for SrVOa, which is used in constraint-LDA calculations, 
structure corresponds to the canonical-bands approximation^- The Fermi level is at zero energy. 
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FIG. 17: A schematic view on the change of the p-d hybridization in the 0(2p) and t2 g bands of the TM oxides associated with 
the repopulation of the Wannier orbitals at the neighboring transition-metal sites: 2d n ^ ± d n ~ 1 +d n+1 . Left panel corresponds 
to the ground-state configuration (2d n ). In the right panel, the removal (addition) of an electron from (to) the Wannier orbital 
in the i2 9 part of the spectrum is simulated by the shift of these orbitals relative to the Fermi level (shown by dot-dashed line). 
Around each transition-metal site, it changes the Coulomb potential, which controls the distribution of the d-states as well as 
the degree of their hybridization between TM(d) and 0(2p) states. Generally, the removal of an electron from (or the addition 
of an electron to) the Wannier orbital is partially compensated by the change of the amount of d-states, which is admixed into 
the 0(2p) band. This transfer of the spectral weight works as a very efficient channel of screening of local Coulomb interactions 
in the transition-metal oxides. 



To some extent, the true screening can be simulated by imposing certain constraints on the form of RPA polarization 
function P H . Namely, one can expect certain cancellation of contributions coming from Wannier orbitals centered at 
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different TM sites at the intermediate (i.e. oxygen) sites (see Fig. I18|) . Therefore, we believe that it is more physical 




FIG. 18: A cartoon picture showing two Wannier orbitals centered at neighboring transition-metal sites (T) and having tails 
spreading to the intermediate oxygen (O) site. The charge disproportionation associated with the reaction 2(d n )^d" +Sn +d"~ s " 
will affect mainly the transition-metal sites. The change of the electron density at the intermediate oxygen site will be largely 
cancelled out. 



to take into account only those contributions to the polarization function which are associated with the TM sites, 
and to suppress contributions associated with intermediate sites. This makes some difference from the conventional 
RPA,2L22iiS which constitutes the basis of the GW method4ii2i We expect this scheme to work well for the on-site 
Coulomb interactions. However, the effect of hybridization on the intersite Coulomb interactions remain an open 
and so far unresolved problem. Some estimates of these effects will be given in Sec. IVII Dl using the c-LDA method. 
The calculations suggest that the intersite interactions are screened very efficiently by the change of hybridization. 
Therefore we speculate that for the considered compounds, the effective Coulomb interactions between different TM 
sites are small and can be neglected. 

The analytical expression for P H can be obtained from Eq. (|27|l by considering the perturbation-theory expansion 
for the wavefunctions with the fixed occupation numbers^ The time-dependent perturbation theory, corresponding 
to the external perturbation SV cx te~ luJt , yields in the first order: 



P. 



H 



(-) = £■ 



[m - n^d^dpid^dsj 



+ Si + iS(rii — rij) ' 



(29) 



where d 1 i=((j) 1 \ipi) is the projection of LDA eigenstate ipi onto one of partial d-waves </> 7 , belonging to the TM site, 
and i and j are the joint index, incorporating spin and band indices as well as the position of k-point in the first 
Brillouin zone. In this notations, the matrix multiplication in Ea. l28l implies the convolution over two indices. Namely, 
the matrix element of the product uP H is given by (uP H ) a 0-yS = Y^ui/ u a p fjLU P^ lS . Note that all transitions in Eq. 
are allowed only between occupied and empty bands. 

In the next Sections, we will discuss different contributions to the screening of Coulomb interactions more in details. 



B. Results of Constraint-LDA Calculations 



The results of conventional constraint-LDA calculations for TM oxides have been widely discussed in the 
literature i^SLS&SSiSi Here we only illustrate the main idea and show some basic results using SrV03 as an exam- 
ple. 

The calculations are performed in the supercell geometry, in which the number of atomic 3d electrons is modulated 
around the "ground-state" value n=l according to the formula: 

fin = n + 5n cos(kR). 

Corresponding values of interaction parameter Uk, calculated in several different points of the Brillouin zone, are listed 
in Table Ixl Then, we map Wk onto the model 
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TABLE X: The values of Mk (in eV) obtained in constraint-LDA calculations for SrVC>3 in different points of the Brillouin zone: 
r=(0,0,0), X=(7r/o) (1,0,0), M=(7r/o)(l,l,0), and R=(7r/a)(l, 1, 1). The right column shows the weights of these k-points 
used in the process of integration over the Brillouin zone. 



k 




w k 


r 







X 


10.6 


3/7 


M 


9.9 


3/7 


R 


9.5 


1/7 



itk = u — v cos(kR), 

R 

and extract parameters of on-site (it) and NN (v) interactions after integration over the Brillouin zone. We note that 
in the present context the T-point result has no physical meaning, as it corresponds to the transfer of an electron 
to the same atomic site. Therefore, we exclude it in the process of integration, and recalculate the weights of other 
k-points using the symmetry arguments. The new weights are shown in Tabled This yields the following parameters: 
m=10.1 eV and v— 1.2 eV. Similar calculations for intra-atomic exchange coupling yield j = 1 . eV. 

For comparison, the values of bare Coulomb and exchange integrals, calculated on the atomic V(3d) wavefunctions 
are 21.7 and 1.2 eV, respectively. Thus, in the c-LDA scheme the effective Coulomb interaction is reduced by factor 
two. The intra-atomic exchange interaction is reduced by 20%. As we will see in the next section, the Coulomb 
interaction will be further reduced by relaxation of the WFs due the change of hybridization. 



C. The Role of Hybridization 

Because of hybridization, the d-states of the TM sites may have a significant weight in other bands. A typical 
situation for the series of TM oxides is shown in Fig. ^ where besides the i 2 g-band, the d-states contribute to the 
TM e g as well as the 0(2p) bands. On both sides of the reaction 2{d n )^±d n+1 +d n ~ l , the WFs constructed for the 
t2 g bands should be orthogonal to other bands. As it was already pointed out in Sec. I VI I Al this mechanism is 
responsible for an additional channel of screening of the on-site Coulomb interaction associated with the change of 
this hybridization. The corresponding contribution can be evaluated using the Dyson equation l|28|l and taking the 
matrix of Coulomb interactions u obtained in c-LDA as the starting interaction. Then, the relevant expression for 
the polarization matrix is given by Eq. I|29() . 

According to the electronic structure of the TM oxides, one can identify three main contributions to the po- 
larization matrix P H , associated with the following inter-band transitions: 0(2p)— >TM(e 9 ), 0(2p)— >TM(i2 g ), and 
TM(i2 9 )— >TM(e 9 ). Meanwhile, all transitions between t2 g bands should be switched off, in order to avoid the double 
counting of these effects in the process of solution of the Hubbard models As it was already pointed out in Sec. lVIII 
this procedure is similar to the setting h=c for the auxiliary WFs. 

Details of static screening, corresponding to uj—0, are explained in Fig.^|for SrVOa. It is convenient to introduce 
three Kanamori parameters^ the intra-orbital Coulomb interaction U=U xy xy xy xy , the inter-orbital interaction 
W=U xy xy yz yz , and the off-diagonal (exchange-type) interaction J=U xy yz xy yz . In addition to the total value of 
U, we calculate intermediate interactions corresponding to each type of transitions in the polarization matrix (and 
neglecting the other two). The screening caused by the change of the hybridization appears to be very efficient. So, 
by going from c-LDA to RPA the intra-orbital interaction U is reduced from 11.2 to 2.5 eV (i.e., by factor four and 
even more). The main contribution to this screening comes from the 0(2p)^V(e g ) and 0(2p)-^V(t2 g ) transitions in 
the polarization functions. In the cubic perovskites, the direct interaction between Y{t2 g ) and V(e s ) bands plays a 
minor role and can be neglected. 

Matrix elements of the (total) polarization function are displayed in Fig. [2U] The largest contribution comes from 
the site-diagonal elements of the type P(oj)=P t ^ aaa (uj)=P l ^ l3 p a (u}). Other contributions are considerably smaller. The 
static polarization P=P(0) is about —0.12 eV -1 . This is the large value because the renormalization of the on-site 
Coulomb interaction in the multi-orbital systems is controlled by the parameter MP, rather than P (see Appendix). 
The prefactor M stands for the total number of orbitals per one TM site (M=3 for t^g systems). Therefore, —uMP can 
be estimates as 3.6, and the situation appears to be close to the strong-coupling regime. Then, the effective interaction 
is not sensitive to the exact value of the parameter u, which is used as the starting point in RPA calculations. For 
example, had we started with the bare Coulomb interaction, which exceed the c-LDA value by factor two and even 
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FIG. 19: Left panel shows the local density of states of SrVC>3 with the notation of the main inter-band transitions which 
contribute to the polarization function in RPA: 0(2p)— >V(e g ) (1), 0(2p)— >V(£2 9 ) (2), and V(t2 9 )— *V(e 9 ) (3). Right panel 
shows the intra-orbital interaction IA in c-LDA (U=U+jJ)~^ the effect of screening corresponding to each type of transitions 
in the polarization function, and the final value of U in RPA which incorporates all three transitions. All RPA results are for 
lu=0. 
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FIG. 20: The diagonal matrix element P xy xy xy xy of the polarization function for SrVOa in the real space: the black curve 
shows the site-diagonal part, the light blue curve shows nearest-neighbor elements in the xy-pl&ne. The real and imaginary 
parts are shown by solid and dashed curves, respectively. 

more, we would have obtained U=2.7 eV, which is close to 2.5 eV derived by starting with c-LDA. This justifies 
some approximations discussed in Sec. IVH Al particularly the use of fast but not extremely accurate c-LDA for some 
channels of screening. 



D. RPA versus Constraint-LDA for t2 g Electrons 

In this section we briefly return to the problem considered in Ref. |3^ and reinterpret some results obtained in that 
work in the light of present RPA approach. The basic idea of Ref. |3^ was to evaluate the effective Coulomb interaction 
for the series of TM perovskite oxides in the framework of c-LDA, which would incorporate the screening by itinerant 
TM(e g ) electrons. Since for the considered type of screening, RPA has many similarities with c-LDA, the present 
section can be also regarded as a test for these two approaches. 

However, it is important to remember that several basic assumptions of Ref.[35| were different from the present work. 
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So, the TM(t2 g ) states in Ref. |35| were totally decoupled from the rest of the electronic states by switching off the 
matrix elements of hybridization in the LMTO method, and only the TM(e g ) states were allowed to hybridize. The 
corresponding electronic structure of SrVC>3 is shown in Fig. [21 In terms of RPA polarization function, this means 
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FIG. 21: Electronic structure of SrVC>3 in canonical-bands approximation for the V(t2 S ) states^- the contributions of V(e 9 ) 
states are shown by arrows. Other notations are the same as in Fig.Q 

that the only allowed contributions to the screening in Ref. l35l were due to the 0(2p)— >V(e ff ) inter-band transitions 

(type 1 in Fig. 033). 

Results of such c-LDA calculations for SrVC-3, which have been performed along the same line as in Sec. IVII Bl 
are summarized in Table [XTI After the Fourier transformation to the real space, we obtain the following parameters 



TABLE XI: Results of constraint-LDA calculations for [7k between ti g states in SrVC>3 obtained after including the e 9 states 
to the screening of t2 g interactions. Other notations are the same as in Table IXl 
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of on-site and NN interactions: U =3.4 eV and V=0.3 eV. This value of U appears to be in a reasonable agreement 
with the final [7=2.5 eV, extracted from RPA (Fig. I19J1 . However, c-LDA employs an additional atomic-spheres 
approximation. Therefore, for a proper comparison with RPA, one should use the same level of approximation and 
suppress all nonspherical interactions in the matrix of Coulomb interactions u, which is used as the starting point in 
RPA. In this approximation, and considering only the 0(2p)— >V(e g ) transitions in the polarization function, we obtain 
C/=3.6 eV, which is close to the c-LDA value obtained using the method proposed in Ref. |3f|. The small difference 
is caused by different approximations used for treating the intersite Coulomb interactions, which were neglected in 
RPA and taken into account in c-LDA. For comparison, the total value of U obtained in RPA after neglecting the 
nonsphericity effects is only 1.6 eV. 

In summarizing this section, there is a reasonable agreement between results of RPA calculations and the c-LDA 
approach proposed in Ref. l35l However, the agreement is somewhat fortuitous because this c-LDA takes into account 
only one part of the total screening, corresponding to the 0(2p)^TM(e g ) transitions in the polarization function. 
The error caused by this approximation is partially compensated by the atomic-spheres approximation supplementing 
the c-LDA scheme. 

The c-LDA calculations give some idea about the effect of hybridization on the screening of intersite Coulomb 
interactions. The screening appears to be very efficient. So, by taking into account only the 0(2p)— >V(e g ) transitions, 
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the NN interaction is reduced from 1.2 to 0.3 eV. We expect this value to be further reduced by including other types 
of transitions in the polarization function. Thus, for the considered compounds, the effective Coulomb interaction 
between different TM sites seems to be small and can be neglected. 



E. Doping-Dependence and Kanamori Rules for Cubic Compounds 



In this section we discuss the effects of electron/hole doping on the static Coulomb interactions in SrV03 using the 
rigid-band approximation. Results of such c-LDA+RPA calculations are shown in Fig. |23 versus the total number of 
electrons in the TM(i 2g ) band, n t2g - We monitor the behavior of three Kanamori parameters^ U, W , and J . 
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FIG. 22: Doping-dependence of the effective Coulomb interactions in SrVOa. 

The Coulomb interactions reveal a monotonic behavior as the function of doping. The screening is the most 
efficient when the t^g band is empty (nj 2 =0). The situation corresponds to SrTiOs. In this case all 0(2p)— >TM(i2 9 ) 
transitions contribute to the screening in RPA (see Fig. I19fl . This channel of screening vanishes when the t 2g band 
becomes occupied {n t2g —6). Then, the only possible screening is associated with the 0(2p)— >TM(e 9 ) transitions in 
the polarization function, and the effective Coulomb interaction becomes large. 

The screening of off-diagonal matrix element J practically does not depend on doping. Therefore, the well know 
Kanamori rule, U=U'+2J', which was originally established for atoms, works well in the cubic compounds, even after 
the screening of ti g interactions by other electrons. 

The present result also supports an old empirical rule suggesting that only the Coulomb integral U is sensitive to 
the crystal environment in solids. The nonspherical interactions, which are responsible for Hund's first and second 
rules, appears to be much closer to their atomic valuesi^il 



F. Frequency-Dependence 

We have shown that the change of hybridization plays a very important role and strongly reduces the static value 
of U. However, this is only one part of the story because the same effect implies the strong frequency-dependence of 
the effective interaction, as it immediately follows from the Kramers-Kronig transformation in RPA4i 

ReUjuj) = u--V f^y^l. (30) 

Indeed, the difference [u— Ret/ (a;)] at u>=0, which for the diagonal matrix elements of SrV03 is about 8.7 eV, should 
be related with the existence of the finite spectral weight of |Imt/(o;)| at finite lo. These dependencies are shown in 
Fig. HO for SrV0 3 . 
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The high-frequency part of U can also contribute to the low-energy part of the spectrum through the self-energy 
effect j22i The latter can be evaluated in the GW approximation^^ where the self-energy is given by the convolution 
of U(uj) with the one-particle Green function G(cu) for the t2 g band: 

t(u) = f dui'G(ui + ui')U(ui'). 
2n J 

For SrV03, the diagonal matrix element of £(w) is also shown in Fig. 1231 The low-frequency part of ImE is small 
and can be neglected, while ReS mainly contributes to the renormalization factor: 

Z aP = [l-dReE^/dul^y 1 . 
The latter is estimated as 0.8, for the diagonal matrix elements. 



G. Lattice Distortion, Formal Valency and Screening 

In Sec. lVIBl we already pointed out that the degree of hybridization between atomic TM(t2 9 ) and 0(2p) states can 
differ substantially for different TM oxides. A typical example is two isoelectronic perovskites: SrV03 and YTi03. 
The TM(t 2g )-0(2p) hybridization is stronger in SrV03, because of two reasons: 

(i) A direct proximity of 0(2p) and V(<2g) bands in SrVC"3, which is expected for tetra-valent compounds; 

(ii) A strong orthorhombic distortion observed in YTi03, which generally deteriorates the Ti(t2 g )-0(2p) hybridization. 
Therefore, it is reasonable to expect a very different screening of on-site Coulomb interactions in these two com- 
pounds. This idea is nicely supported by results of RPA calculations shown in Fig. [2] The static U is larger YTi03. 
For example, the diagonal matrix element of U is about 3.4 eV, against 2.5 eV in SrV03. This is despite the fact 
that c-LDA has an opposite tendency. Generally, the value of u in c-LDA is expected to be larger for SrV03 rather 
than for YTi03, due to different shape of the atomic 3d-wavefunctions in the three- and tetra-valent compounds^ 
So, c-LDA yields u— 8.9 and 10.1 eV, correspondingly for YTi03 and SrV03. These parameters have been used as 
the starting point in RPA, which results in an opposite trend for the static U. 

Therefore, the RPA screening is more efficient in SrV03, which is consistent with stronger TM(t2 g )-0(2p) hybridiza- 
tion in this compound. On the other hand, the frequency-dependence of U is weaker in YTi03, as it immediately 
follows from the Kramers-Kronig transformation (|30l) . 

Finally, because of different hybridization of the t2 g orbitals in YTi03, the diagonal matrix elements of U are also 
different (see Fig. I24f) . In this case, there is some deviation from the Kanamori rules. The difference is small (about 
0.07 eV for diagonal matrix elements of U at w=0). However, it may play some role in more delicate applications, 
such as the orbital magnetism in solids, for example. 34 
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FIG. 24: Diagonal matrix elements of screened Coulomb interactions for SrV03 (left) and YTiOa (right). Three different lines 
in the case of YTi03 show the behavior of three diagonal matrix elements in the local coordinate frame given by Eq. [1201 . 
The insets show the high-frequency part of U. The lines corresponding to different ti g orbitals in the case of YTiOa are 
indistinguishably close. 



H. Two Models for V 2 3 

In Sec. I VI iJl we introduced two possible models for the kinetic-energy part of V2O3: the "five-orbital model" treats 
all V(3d) bands on an equal footing, while the "three-orbital model" is limited by twelve Vfeg) bands, close to the 
Fermi level. We argued that parameters of the kinetic energy can be different for these two models. The same is true 
for the effective Coulomb interactions. The RPA provides a very transparent explanation for this difference, which is 
based on the following arguments. Recall that RPA incorporates the screening of on-site Coulomb interactions caused 
by relaxation of the wavefunctions. This relaxation is treated analytically, using regular perturbation theory for the 
wavef unctions, which results in Eq. (|29[1 for the polarization function^ 

Since the V(e g ) band is eliminated in the three-orbital model, the effective U should include the screening caused 
by that change of the wavefunctions, which is formulated in terms of transitions between V(t2 g ) and V(e g ) bands in 
the perturbation-theory expansion. In the five-orbital model, the V(e g ) band is included explicitly. Therefore, the 
relaxation caused by possible interactions between V(i2g) and V(e g ) bands will be automatically taken into account 
in the process of solution of the five-orbital model, and we should get rid of this parasitic screening at the stage 
the construction of the model Hamiltonian. Thus, the tree-orbital model will include an additional screening of on- 
site interactions, caused by the V(t2 S )— »V(e 9 ) transitions in the polarization function, which does not appear in the 
five-orbital model. 

This screening appears to be rather efficient (unlike in cubic perovskites considered in Sec lVII C|) . So, the diagonal 
matrix element of static t2 g interactions is about 3.2 and 3.9 eV, for the three- and five-orbital model, respectively 
(Fig. I25|l . On the other hand, according to the Kramers-Kronig transformation l|30|) . the frequency-dependence of 
the effective interaction is more important in the three-orbital model. Because of different hybridization of tig and 
e g states in the five-orbital model, the e g interactions are screened more efficiently (the static interaction between e g 
orbitals is about 3.7 eV). 



VIII. SUMMARY, OPEN QUESTIONS, AND COMPARISON WITH OTHER METHODS 

The ultimate goal of this work was to make a bridge between first-principle electronic structure calculations and the 
universe of Hubbard parameters for strongly-correlated systems. We have presented a comprehensive analysis of the 
problem, by starting with the brief description of the ASA-LMTO method for electronic structure calculations and 
ending up with realistic parameters of the kinetic-energy and the Coulomb interactions for the series of TM oxides 
obtained on the basis of this LMTO method. A particular attention has been paid to the analysis of microscopic 
processes responsible for the screening of on-site Coulomb interactions in oxide compounds. 

Our strategy consists of three steps: 
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FIG. 25: Diagonal matrix elements of screened Coulomb interactions obtained in the three- (left) and five-(right) orbitals 
models for V2O3. The notations of orbitals are given in Sec. IVI Dl Other notations are the same as in Fig. 1241 



(i) Derivation of the kinetic-energy part of the Hubbard model from the single-particle electronic structure in LDA, 
using the downfolding method. We have also considered corrections to the crystal-field splitting caused by nonspheric- 
ity of electron-ion interactions, beyond the conventional atomic-spheres-approximation. 

(ii) Construction of the Wannier functions using results of the downfolding method. At this stage we closely follow the 
idea of LMTO method, and construct the WFs as the LMTO basis functions, which after applying to the Kohn-Sham 
Hamiltonian in the real space generate the matrix elements of the kinetic energy obtained in the downfolding method. 

(iii) Calculation of screened Coulomb interactions using the concept of auxiliary WFs. The latter are defined as the 
Wannier orbitals for which the kinetic-energy part is set to be zero. This construction allows to avoid the double count- 
ing of the kinetic-energy term, which is included explicitly in the Hubbard model. The screened Coulomb interactions 
are calculated on the basis of a hybrid approach, combining the conventional constraint-LDA with the random-phase 
approximation for treating the hybridization effects between atomic TM(3g?) and 0(2p) orbitals. The latter play a 
very important role and yields a strong renormalization of the effective Coulomb interaction for isolated ti g band. It 
also explains a strong material-dependence of this interaction, which is sensitive to the crystal environment in solids, 
the number of t2 g electrons, the valent state of the TM ions, etc. 

Taking into account a wide interest to the construction of effective lattice fermion models from the first principles, 
we would like to make a brief comparison with other works on a similar subject. 

Majority of methods start with the construction of the WFs, which are then used as the basis for calculations of 
the parameters of the kinetic energy and the Coulomb interactions. This is different from our approach, where we 
start with the kinetic-energy part, and only after that construct the WFs for a given set of parameters of the kinetic 
energy. We believe that such an order is extremely important, as it allows us to control the contributions of the 
kinetic energy to the WFs and the Coulomb interactions. 

Among recent works, a considerable attention is paid to the method of Marzari and Vanderbilt, because it allows 
to control the spacial extension of the WFs. Very recently, Schnell et al. applied this method to calculations of the 
parameters of the Hubbard Hamiltonian for the series of 3d transition metals4& In each k-point of the Brillouin zone, 
they constructed the WFs from all 16 bands of the LDA Hamiltonian, corresponding to the 4s4p3<i4/ LMTO basis. 
Therefore, the total number of Wannier orbitals was also 16. The coefficients of this expansion has been chosen so to 
minimize the square of the position operator, (r 2 )=(W\r 2 \W) . The obtained WFs were indeed well localized, and the 
parameter U estimated for the 3c? bands was very close to the atomic value (about 25 eV). However, the corresponding 
Wannier basis set is too large, that does not make a big difference from the original LMTO basis set, for which one can 
also introduce a localized (tight-binding) representation. 1 - From the view point of numerical solution of the Hubbard 
model using modern many-body techniques, it is still hardly feasible to work in the basis of 16 Wannier orbitals per 
one TM site, while the simplest Hartree-Fock approximation is definitely not sufficient for the transition metals4& It 
would be interesting to see how this method will work for the TM oxides, considered in the present work, where the 
physical basis set is limited by three Wannier orbitals per one TM site. For example, is it possible to construct the 
localized Wannier orbitals for isolated t^g bands in the TM oxides, which would be as good as the Wannier orbitals 
derived for all bands? The problem is that when the number of bands decreases, the number of variational parameters 
for the optimization of the Wannier orbitals will also decreases. Therefore, the Wannier orbitals will generally become 
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less localized, as it was demonstrated in our work for V2O3. Also, when the problem is formulated in a reduced Hilbert 
space of bands closest to the Fermi level, it is very important to consider the screening of Coulomb interactions by 
other bands, which comes from relaxation of the wavefunctions. These effects are beyond the scopes of the work of 
Schnell et al., who only considered the bare Coulomb interactions. 

In order to construct the WFs, Ku et al. employed the projector-operator scheme^ which is basically the initial 
step of the method of Marzari and Vanderbilt, prior the optimization. 13 In this scheme, each Wannier orbital is 
generated by projecting a trial wavefunction, \g), onto a chosen subset of bands (in our case, tig bands): 

\W) = \A)(A\g)- 

Since such orbitals are not orthonormal, the procedure is followed by the numerical orthonormalization, similar to the 
one described in Sec. [H] These WFs have been used as the basis for the construction of the low-energy Hamiltonian 
for the series of cuprates, like La4f3a2Cu20io, where the Wannier basis consisted of a single orbital centered around 
each Cu site. A weak point of this approach is that it is difficult to assess the spacial extension of the WF, which 
crucially depends on the trial wavefunction, and can be affected by the orthonormalization. In some sense, the result 
strongly depends on authors' intuition on how they choose the trial wavefunction. For example, the bare Coulomb 
interaction obtained in Ref. \^ in the basis of their WFs was only 7.5 eV, which is much smaller than the atomic 
Coulomb integral for the Cu(3d) orbitals. This means that the WFs are not well localized. It is not clear at present, 
whether this is a result of the bad choice of the trial wavefunction, or there is a more fundamental problem related 
with the fact that a more compact representation for the WFs simply may not exist in this case. Note, that apart 
from the Berry phase, there is no further parameters available for the optimization of the WFs in the single-orbital 
case. 

The derealization of the WFs gives rise to appreciable direct exchange interactions operating between different Cu 
sites^?. This result, however, rises additional questions. Note, that the kinetic part of the Hubbard model and the 
WFs are evaluated in LDA, where the exchange-correlation potential is set to be local. In principle, non-local effects 
can be already incorporated in LDA, through the renormalization of parameters of the kinetic energy and the local 
interactions^ Therefore, it is not clear whether the nonlocal exchange interactions should be regarded as independent 
parameters of the Hubbard Hamiltonian or not. 

Finally, Ku et al. calculated only bare Coulomb interactions. They did not consider the screening of these interac- 
tions caused by relaxation effects, which are extremely important. 

Anisimov et. al. employed a similar approach for the analysis of spectroscopic properties of TM oxides^i They 
extracted only the kinetic-energy part of the Hubbard Hamiltonian, using the WFs constructed in the LMTO basis, 
and treated the Coulomb interaction U as a parameter. 

A completely different strategy has been proposed by Andersen et. al, on the basis of their order- N muffin-tin 
orbital (NMTO) method, which is an extension of the LMTO method^ From the very beginning, they construct the 
NMTO basis functions in certain energy interval as the WFs of the original KS Hamiltonian. Pavarini et al. applied 
this method to the series of al 1 perovskitesS At present, it is not clear how the nonuniqueness of the WFs is reflected 
in the construction of the NMTO basis set. Obviously, such a basis set is also not unique, and there is some freedom 
left for the localization of the Wannier orbitals, which does not seem to be well controlled. Generally, our transfer 
integrals for the d 1 perovskites seems to be more localized and our crystal-field splitting is smaller. For example, had 
we relaxed the constraint condition for the construction of the "heads" of the WFs, based on the diagonalization of 
the density matrix (|I9|I . our conclusion would have been also different: the crystal-field splitting would increase, but 
the transfer integrals would become less localized. 

Pavarini et al. did not calculate the Coulomb interactions. Instead, they used C/~5 eV as a parameter, with the 
reference to the photoemission data^ However, the photoemission data are typically interpreted in the cluster model, 
which treats the 0(2p) band explicitly. For the isolated t2 g band, the effective interaction should include an additional 
renormalization coming from the relaxation of the 0(2p) band, which is eliminated in the £2g-model. Therefore, the 
value of the effective U should be smaller. 

Finally, due to unknown for us reason, there is a substantial difference of the parameters of ti g bandwidth (Wt 2g ) 
between our work and Ref. l53l even for cubic SrV03. The parameters reported by Pavarini et al. are generally 
overestimated by about 30%S Our conclusions about magnetic properties of YTi03 and LaTi03 are also differentiii 
Details will be presented in a separate paper. 
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APPENDIX: ORBITAL DEGENERACY AND SCREENING IN RPA 



In this Appendix we consider the screening of on-site Coulomb interactions in RPA for an M-orbital system. All 
orbitals are supposed to be equivalent. For simplicity we neglect small nonsphericity of bare Coulomb interactions. 
Then, the nonvanishing matrix elements of the Coulomb interactions are u aa pp = u, where a (/3)=1, . . . ,M. They 
can be presented in the form u=ul, where I is the MxM matrix, consisting of only the units: 



I = 



f 1 1 
1 1 



1\ 



Vi iJ 



The part of the polarization polarization matrix (|29[> . which can interact with the matrix u, is assumed to be diagonal: 
Paapp = P5 a f3- The assumption is justified for cubic perovskites, where different ti g orbitals belong to different bands. 
For other compounds it can be regarded as an approximation, which does not change our qualitative conclusion. Hence, 
for the screened Coulomb interaction l|28[) we have: 



U=[l- uP^u = ^{uP) n i n uI. 



(A.l) 



Since I n+1 =M n I, Eq. ljA~T|) can be converted to 



U = ^2(uMP) n ui = - 



n=0 



-uMP 



-I. 



This means that in the multi-orbital systems, the renormalization of the Coulomb repulsion is more efficient as it is 
controlled by the quantity (MP), where the prefactor M stands for the number of orbitals. 

In the strong coupling limit, -«MP>1, the effective interaction is U=—(MP)~ 1 I, which does not depend on the 
value of bare Coulomb interaction u. 
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